State-Space Model Analysis

Author

Tzu-Yao Lin

Published

February 21, 2024

1 Dataset

This dataset comes from Koval and his colleagues’ (2013) study, in which they use a novel paradigm and analytic approach to study the dynamic process of depression. In past studies, depression is associated with higher average levels of negative affect (NA) and lower average levels of positive affect (PA). However, few studies untangle the relationship between the pattern of affective fluctuations and the dynamics of depression. Therefore, in their study, they collected 99 (the final number) subjects’ depression symptoms by the Center for Epidemiologic Studies Depression Scale and (CES-D) and daily experiences of affect ratings 10 times/day for 7 days using the experience sampling method (ESM).

Affect ratings were made on:

  • PA = mean(happy, relaxed), and

  • NA = mean(sad, depressed, anxious, angry).

The scale of each affect rating is from 1 (not at all) to 100 (very much).

library(tidyverse)
theme_set(theme_bw(base_size = 14))
library(lubridate)

rawdata <- read_tsv("data_1beep_no1st beep_annette.tsv")
rmarkdown::paged_table(rawdata)
data <- rawdata %>% 
  mutate(Pos = PA, 
         Neg = `NA`,
         Date_Time = ymd_hms(str_glue("{year}-{month}-{day} {hour}:{min}:{sec}"), tz = "CET"),
         Date = as_date(Date_Time),
         Time = hms::as_hms(Date_Time),
         WDay = wday(Date, label = TRUE),
         Subject = factor(cumsum(PpID != lag(PpID, default = 0))), 
         .keep = "none") %>% 
  group_by(Subject) %>% 
  mutate(Day = factor(cumsum(Date != lag(Date, default = origin)))) %>% 
  group_by(Subject, Date) %>% 
  mutate(Moment = factor(1:n())) %>% 
  ungroup() %>% 
  pivot_longer(cols = c(Pos, Neg), names_to = "Affect", values_to = "Score") 

rmarkdown::paged_table(data)

Nevertheless, I found some discrepancies between the data I got and the descriptions in the literature.

  1. The number of subjects in my data is 97;
  2. The sampling time points are not always 7 (days) and 10 (times). There are many missing values.

All analyses were conducted using the maximum available sample size: for analyses involving ESM data, n=95; for all other analyses, n=99.

t1 <- data %>% 
  group_by(Subject, Affect) %>% 
  summarise(num_day = length(unique(Day)),
            num_obser = n(),
            num_missing = sum(is.na(Score)),
            num_obser_valid = num_obser - num_missing) %>% 
  ungroup()
rmarkdown::paged_table(t1)
summary(t1)
    Subject       Affect             num_day        num_obser    
 1      :  2   Length:194         Min.   :7.000   Min.   :62.00  
 2      :  2   Class :character   1st Qu.:7.000   1st Qu.:66.00  
 3      :  2   Mode  :character   Median :7.000   Median :66.00  
 4      :  2                      Mean   :7.031   Mean   :66.75  
 5      :  2                      3rd Qu.:7.000   3rd Qu.:67.00  
 6      :  2                      Max.   :8.000   Max.   :78.00  
 (Other):182                                                     
  num_missing     num_obser_valid
 Min.   : 0.000   Min.   :39.00  
 1st Qu.: 3.000   1st Qu.:58.00  
 Median : 6.000   Median :61.00  
 Mean   : 6.737   Mean   :60.02  
 3rd Qu.: 9.000   3rd Qu.:64.00  
 Max.   :29.000   Max.   :73.00  
                                 
table(t1$num_day)

  7   8 
188   6 
t2 <- data %>% 
  group_by(Subject, Affect, Day) %>% 
  summarise(num_moment = n(), 
            num_missing = sum(is.na(Score)), 
            num_valid = num_moment - num_missing) %>% 
  ungroup()
rmarkdown::paged_table(t2)
summary(t2)
    Subject        Affect               Day        num_moment    
 52     :  16   Length:1364        1      :194   Min.   : 5.000  
 70     :  16   Class :character   2      :194   1st Qu.: 9.000  
 87     :  16   Mode  :character   3      :194   Median :10.000  
 1      :  14                      4      :194   Mean   : 9.494  
 2      :  14                      5      :194   3rd Qu.:10.000  
 3      :  14                      6      :194   Max.   :20.000  
 (Other):1274                      (Other):200                   
  num_missing        num_valid     
 Min.   : 0.0000   Min.   : 0.000  
 1st Qu.: 0.0000   1st Qu.: 8.000  
 Median : 1.0000   Median : 9.000  
 Mean   : 0.9582   Mean   : 8.536  
 3rd Qu.: 1.0000   3rd Qu.:10.000  
 Max.   :10.0000   Max.   :11.000  
                                   
table(t2$num_moment)

  5   6   7   8   9  10  11  14  20 
 10  20 138  42 188 864  98   2   2 
table(t2$num_valid)

  0   2   3   4   5   6   7   8   9  10  11 
  2   4   2  30  46  89 147 208 348 448  40 

2 Data Illustration

library(tsibble)
subdata <- data %>% 
  filter(Subject %in% c(9, 20, 36, 57, 76, 85)) %>% 
  mutate(DateTime = as_datetime(ymd_hms(paste(as_date(as.double(Day)), 
                                              as.character(Time))))) %>% 
  as_tsibble(key = c(Subject, Affect), 
             index = DateTime) 

pos_neg_color <- rev(scales::hue_pal()(2))

ggplot(subdata, aes(x = DateTime, y = Score, color = Affect)) + 
  geom_line() + geom_point() +
  scale_y_continuous(limits = c(0, 100)) +
  scale_x_datetime(breaks = as_datetime(1:7 * 86400),
                   labels = paste("Day", 1:7),
                   limits = as_datetime(c(1, 8) * 86400)) +
  scale_color_manual(values = pos_neg_color) +
  facet_grid(rows = "Subject") 
Figure 1: Illustration of the Positive and Negative Affect scores for different time points for a subgroup of subjects

Some observation:

  • The positive affect scores generally are larger than the negative affect scores. Both scores are negatively correlated. But there are some exception cases.

  • The sampling time points are quite different between subjects.

  • The data was collected majorly in the latter half of the day.

  • The sampling intervals are not fixed within days and subjects.

  • There are a lot of missing data.

3 State-Space Model

3.1 Single Level Model

Here, I only fitted one subject specifically.

3.1.1 Model Specification

Here, I follow Schuurman & Hamaker (2019) Measurement Error Vector Autoregressive with Order 1 Model (MEVAR(1)) without a multilevel setting. It means I will fit each subject separately.

  • Measurement equation

\[\begin{pmatrix}y_{1it} \\ y_{2it}\end{pmatrix} = \begin{pmatrix}\mu_{1i} \\ \mu_{2i}\end{pmatrix} + \begin{pmatrix}\theta_{1it} \\ \theta_{2it}\end{pmatrix} + \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \]

\[ \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \sim \mathcal{N} \left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \mathbf{R}_i = \begin{pmatrix} \sigma^2_{\epsilon 11i} & \sigma_{\epsilon 12i} \\ \sigma_{\epsilon 12i} & \sigma^2_{\epsilon 22i}\end{pmatrix}\right) \]

  • State equation (state space model representation)

\[ \begin{pmatrix}\theta_{1it} \\ \theta_{2it}\end{pmatrix} = \begin{pmatrix} \phi_{11i} & \phi_{12i} \\ \phi_{12i} & \phi_{22i} \end{pmatrix} \begin{pmatrix}\theta_{1it-1} \\ \theta_{2it-1}\end{pmatrix} + \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \]

\[ \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \sim \mathcal{N} \left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \mathbf{Q}_i = \begin{pmatrix} \sigma^2_{\omega 11i} & \sigma_{\omega 12i} \\ \sigma_{\omega 12i} & \sigma^2_{\omega 22i}\end{pmatrix}\right) \]

where

  • subject index \(i = 1, \dots, N\) and occasion index: \(t = 1, \dots T\)

  • observation vector: \(\mathbf{y}_{it} = \begin{pmatrix} y_{1it} \\ y_{2it} \end{pmatrix}\)

  • latent state (within-person fluctuation): \(\boldsymbol{\theta}_{it} = \begin{pmatrix} \theta_{1it} \\ \theta_{2it}\end{pmatrix}\)

    • initial latent state: \(\boldsymbol{\theta}_{i0} = \begin{pmatrix} \theta_{1i0} \\ \theta_{2i0} \end{pmatrix} \sim \mathcal{N} \left( \mathbf{m}_0 = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \mathbf{C}_0 = \begin{pmatrix} 10^3 & 0 \\ 0 & 10^3 \end{pmatrix} \right)\)
  • person-specific means (trait): \(\boldsymbol{\mu}_i = \begin{pmatrix}\mu_{i1} \\ \mu_{2i}\end{pmatrix}\)

  • measurement errors: \(\boldsymbol{\epsilon}_i = \begin{pmatrix} \epsilon_{1i} \\ \epsilon_{2i} \end{pmatrix}\)

  • autoregressive coefficients: \(\boldsymbol{\Phi}_i = \begin{pmatrix} \phi_{11i} & \phi_{12i} \\ \phi_{12i} & \phi_{22i} \end{pmatrix}\)

  • innovation: \(\boldsymbol{\omega}_i = \begin{pmatrix} \omega_{1i} \\ \omega_{2i} \end{pmatrix}\)

3.1.2 Properties

  • A linear mixed-effect model expression:

We can rewrite the SSM as the LMM. If we combine the measurement equation and the state equation and extend the current time \(t\) to the origin, we have \[ \mathbf{y}_{it} = \boldsymbol{\mu}_{i} + \boldsymbol{\Phi}_i^t\boldsymbol{\theta_{i0}} + \sum_{k=0}^{t-1} \boldsymbol{\Phi}_i^k \mathbf{w}_{i,t-k} + \boldsymbol{\epsilon}_{it} \]

Let \(\mathbf{W}_{it} = \boldsymbol{\Phi}_i^t\boldsymbol{\theta_{i0}} + \sum_{k=0}^{t-1} \boldsymbol{\Phi}_i^k \mathbf{w}_{i,t-k}\), which has a mean \(0\) and a variance \(\mathbf{T}_i\) (see below). So,

\[ \mathbf{y}_{it} = \boldsymbol{\mu}_{i} + \mathbf{W}_{it} + \boldsymbol{\epsilon}_{it} \]

where \(\boldsymbol{\mu}_i\) is a fixed effect, \(\mathbf{W}_{it} \sim \mathcal{N}(\mathbf{0}, \mathbf{T}_i)\) is a Gaussian stochastic process???, and \(\epsilon_{it} \sim \mathcal{N}(0, \mathbf{R})\) is a measurement error. However, there is no random effect in this model under this case.

  • The expectation of \(\mathbf{y}_{it}\)

\[\begin{align*} E[\mathbf{y}_{it}] &= E[\boldsymbol{\mu}_{i}] + E[\boldsymbol{\Phi}_i^t\boldsymbol{\theta_{i0}}] + E[\sum_{k=0}^{t-1} \boldsymbol{\Phi}_i^k \mathbf{w}_{i,t-k}] + E[\boldsymbol{\epsilon}_{it}] \\ &= \boldsymbol{\mu}_{i} + \boldsymbol{\Phi}_i^t \mathbf{m}_0 + 0 + 0 = \boldsymbol{\mu}_{i} \end{align*}\]

  • The variance of \(\mathbf{y}_{it}\)

\[\begin{equation} Var(\mathbf{y}_{it}) = Var(\boldsymbol{\mu}_{i}) + Var(\boldsymbol{\theta}_{it}) + Var(\boldsymbol{\epsilon}_{it}) = 0 + Var(\boldsymbol{\theta}_{it}) + \mathbf{R}_i \end{equation}\]

\[\begin{align} Var(\boldsymbol{\theta}_{it}) = Var(\boldsymbol{\Phi}_i\mathbf{\theta}_{it-1} + \boldsymbol{w}_{it}) = \boldsymbol{\Phi}_i Var(\mathbf{\theta}_{it-1}) \boldsymbol{\Phi}_i^\top + \mathbf{Q}_i \end{align}\]

Under the stationary assumption, we let \(\mathbf{T}_i = Var(\boldsymbol{\theta}_{it})\) for all \(t\).

\[\begin{align*} &\mathbf{T}_i = \boldsymbol{\Phi}_i \mathbf{T}_i \boldsymbol{\Phi}_i^\top + \mathbf{Q}_i \\ \Rightarrow &vec(\mathbf{T}_i) = (\boldsymbol{\Phi}_i \otimes \boldsymbol{\Phi}_i) vec(\mathbf{T}_i) + vec(\mathbf{Q}_i) \\ \Rightarrow &vec(\mathbf{T}_i) = (\mathbf{I} - \boldsymbol{\Phi}_i \otimes \boldsymbol{\Phi}_i)^{-1} vec(\mathbf{Q}_i) \\ \Rightarrow &\mathbf{T}_i = mat((\mathbf{I} - \boldsymbol{\Phi}_i \otimes \boldsymbol{\Phi}_i)^{-1} vec(\mathbf{Q}_i)) \end{align*}\]

Thus,

\[ Var(\mathbf{y}_{it}) = \mathbf{T}_i + \mathbf{R}_i \]

  • The covariance of \(\mathbf{y}_{it}\) and \(\mathbf{y}_{it'}\) (\(t \neq t'\)):

\[\begin{align} Cov(\mathbf{y}_{it}, \mathbf{y}_{it'}) &= Cov(\mu_i + \theta_{it} + \epsilon_{it}, \mu_i + \theta_{it'} + \epsilon_{it'}) \\ &= ????? \end{align}\]

3.1.3 Reliability (subject specific)

  • Based on Schuurman & Hamaker (2019)

\[ R_k = \frac{(\mathbf{T}_i)_{kk}}{(\mathbf{T}_i)_{kk} + (\mathbf{R}_i)_{kk}} \]

where \(k = 1, 2\) for positive and negative affects, respectively.

  • Based on Vangeneugden et al. (2005) and Molenberghs et al. (2007),

There is a problem of representing \(Y\). In Molenberghs et al. (2007), \(\mathbf{Y}_i = (y_{i1}, \dots, y_{iT})^\top\) is a vector collecting observations for all time points. However, in our MEVAR(1) model, \(\mathbf{y}_{it} = (y_{1it}, y_{2it})^\top\) is a vector collecting the positive and negative affects at \(t\). Maybe, I should represent

\[ \mathbf{Y}_i = \begin{pmatrix} y_{1i1} & y_{2i1} \\ \vdots & \vdots \\ y_{1it} & y_{2it} \\ \vdots & \vdots \\ y_{1iT} & y_{2iT} \end{pmatrix} = \mu + ... \]

3.1.4 Bayesian Estimation by Stan

3.1.4.1 Case 1: Fit only for the positive affect

library(cmdstanr)
register_knitr_engine(override = TRUE)
library(posterior)
library(bayesplot)
color_scheme_set("red")
bayesplot_theme_set(theme_bw(base_size = 14))
# library(loo)

The corresponding Stan codes for the one affect of one subject:

ssm0.stan
data {
  int<lower=1> T;
  vector[T] y;
  real m_0;
  cov_matrix[1] C_0;
}

parameters {
  real mu;
  real theta_0;
  vector[T] theta;
  real Phi;
  cov_matrix[1] R;
  cov_matrix[1] Q;
}

model {
  theta_0 ~ normal(m_0, sqrt(C_0[1, 1]));
  
  theta[1] ~ normal(Phi * theta_0, sqrt(Q[1, 1]));
  y[1] ~ normal(mu + theta[1], sqrt(R[1, 1]));
  
  for (t in 2:T) {
    theta[t] ~ normal(Phi * theta[t-1], sqrt(Q[1, 1]));
    y[t] ~ normal(mu + theta[t], sqrt(R[1, 1]));
  }
}

generated quantities {
  vector[T] y_hat;
  
  y_hat = mu + theta;
}

MCMC results:

s20_pos_data <- data %>% 
  filter(Subject == 20, Affect == "Pos") %>% 
  drop_na(Score) %>% 
  pull(Score)

ssm0_data <- list(`T` = length(s20_pos_data),
                  y = s20_pos_data, 
                  m_0 = 0, 
                  C_0 = matrix(10^3))

ssm0 <- cmdstan_model("stan/ssm0.stan")

ssm0_fit <- ssm0$sample(data = ssm0_data, 
                        seed = 20240207,
                        chains = 4,
                        parallel_chains = 4,
                        refresh = 500)

ssm0_fit$save_object(file = "stan/ssm0-fit.RDS")
ssm0_fit <- readRDS("stan/ssm0-fit.RDS")
ssm0_fit$cmdstan_summary() 
Inference for Stan model: ssm0_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.

Warmup took (0.81, 0.85, 0.75, 0.53) seconds, 2.9 seconds total
Sampling took (0.49, 0.50, 0.69, 0.35) seconds, 2.0 seconds total

                 Mean     MCSE  StdDev        5%    50%       95%    N_Eff  N_Eff/s    R_hat

lp__             -370      7.8      29  -4.0e+02   -379  -3.1e+02       13      6.6      1.3
accept_stat__    0.83  4.9e-02    0.25      0.13   0.93      1.00  2.6e+01  1.3e+01  1.1e+00
stepsize__       0.18  5.9e-02   0.083     0.070   0.19      0.30  2.0e+00  9.9e-01  1.5e+14
treedepth__       4.3  3.2e-01    0.85       3.0    4.0       6.0  6.9e+00  3.4e+00  1.3e+00
n_leapfrog__       25  6.6e+00      17       7.0     15        63  6.3e+00  3.1e+00  1.4e+00
divergent__     0.027      nan    0.16      0.00   0.00      0.00      nan      nan      nan
energy__          405  7.9e+00      29       341    414       438  1.4e+01  6.8e+00  1.3e+00

mu                 55     0.13     3.3   5.0e+01     55   6.0e+01      643      318      1.0
theta_0         0.077      1.8      29  -4.8e+01   0.58   4.4e+01      244      121      1.0
theta[1]         0.79     0.22     9.7  -1.6e+01    1.1   1.6e+01     1980      978      1.0
theta[2]         -4.9     0.46     9.1  -2.0e+01   -5.6   1.1e+01      398      196      1.0
theta[3]          2.7     0.24     8.5  -1.2e+01    3.1   1.6e+01     1254      619      1.0
theta[4]          4.3     0.23     8.5  -1.0e+01    4.4   1.8e+01     1407      695      1.0
theta[5]          6.8     0.36     8.9  -8.2e+00    7.3   2.2e+01      621      307      1.0
theta[6]           11      1.1     9.9  -5.7e+00     12   2.6e+01       77       38      1.1
theta[7]          4.8     0.28     8.6  -9.2e+00    5.0   1.9e+01      907      448      1.0
theta[8]          4.8      1.0     9.4  -1.1e+01    5.6   1.9e+01       85       42      1.1
theta[9]          -16      1.5      11  -3.2e+01    -16   2.7e+00       55       27      1.1
theta[10]         -15     0.81      10  -3.1e+01    -16   2.3e+00      163       80      1.0
theta[11]         -19      1.5      11  -3.5e+01    -20   4.8e-01       60       30      1.1
theta[12]         -15     0.85      10  -3.1e+01    -16   2.0e+00      143       71      1.0
theta[13]         -13      1.0      10  -2.9e+01    -14   3.8e+00       96       48      1.1
theta[14]        -1.1     0.20     8.3  -1.5e+01  -0.76   1.2e+01     1680      830      1.0
theta[15]         8.3      1.5      11  -1.0e+01    8.8   2.5e+01       54       27      1.1
theta[16]         -11     0.70     9.3  -2.5e+01    -11   4.8e+00      177       87      1.0
theta[17]         -21      2.0      13  -3.9e+01    -22  -1.1e-01       41       20      1.1
theta[18]         -19      1.9      12  -3.7e+01    -20   9.3e-01       44       22      1.1
theta[19]        0.14     0.53     8.9  -1.6e+01   0.96   1.4e+01      289      143      1.0
theta[20]        -3.6     0.48     8.8  -1.7e+01   -4.2   1.2e+01      344      170      1.0
theta[21]         9.6      1.1     9.6  -6.3e+00     10   2.5e+01       79       39      1.1
theta[22]          10     0.99     9.7  -6.0e+00     11   2.5e+01       96       47      1.1
theta[23]        0.24     0.30     8.7  -1.3e+01  -0.38   1.5e+01      854      422      1.0
theta[24]         9.3      1.0     9.7  -6.8e+00    9.8   2.4e+01       89       44      1.1
theta[25]         4.1     0.60     9.1  -1.1e+01    4.7   1.8e+01      230      113      1.0
theta[26]         -11      1.1      10  -2.6e+01    -12   5.4e+00       82       40      1.1
theta[27]        -8.2     0.27     8.6  -2.3e+01   -8.4   6.0e+00     1000      494      1.0
theta[28]         -16      1.4      11  -3.3e+01    -16   2.5e+00       62       31      1.1
theta[29]        -7.9     0.33     8.8  -2.2e+01   -8.2   6.7e+00      734      363      1.0
theta[30]        0.19     0.32     8.5  -1.5e+01   0.77   1.4e+01      716      353      1.0
theta[31]        -3.3     0.29     8.7  -1.7e+01   -3.8   1.2e+01      910      449      1.0
theta[32]         8.2      1.2     9.7  -7.7e+00    8.7   2.3e+01       69       34      1.1
theta[33]       -0.69     0.22     8.3  -1.4e+01   -1.1   1.4e+01     1432      707      1.0
theta[34]         2.9     0.19     8.6  -1.1e+01    2.7   1.7e+01     1998      987      1.0
theta[35]          21      2.9      16  -3.4e+00     20   4.4e+01       28       14      1.2
theta[36]        -8.4     0.88     9.4  -2.2e+01   -9.1   7.8e+00      114       56      1.0
theta[37]         -16      1.2      11  -3.2e+01    -17   2.0e+00       79       39      1.1
theta[38]         -23      2.2      14  -4.3e+01    -23  -7.2e-02       38       19      1.1
theta[39]        -7.9     0.26     8.7  -2.3e+01   -7.9   6.2e+00     1148      567      1.0
theta[40]        0.97     0.73     9.1  -1.5e+01    1.6   1.5e+01      156       77      1.0
theta[41]         -13      1.4      11  -2.9e+01    -14   4.4e+00       56       28      1.1
theta[42]        -7.4     0.40     8.9  -2.1e+01   -7.9   7.6e+00      491      243      1.0
theta[43]        -1.8     0.19     8.5  -1.6e+01   -2.1   1.2e+01     2062     1018      1.0
theta[44]          10     0.86     9.5  -5.3e+00     11   2.5e+01      122       60      1.0
theta[45]          19      1.6      12  -4.4e-01     20   3.6e+01       51       25      1.1
theta[46]          14     0.79     9.8  -2.1e+00     15   2.9e+01      155       77      1.0
theta[47]          17      1.2      10  -1.8e+00     17   3.2e+01       82       41      1.1
theta[48]          14     0.86     9.9  -2.6e+00     15   3.0e+01      134       66      1.0
theta[49]          16      1.4      11  -2.3e+00     17   3.2e+01       57       28      1.1
theta[50]         8.0     0.71     9.0  -6.8e+00    8.4   2.1e+01      161       79      1.0
theta[51]         -11      1.7      11  -2.8e+01    -11   7.5e+00       43       21      1.1
theta[52]       -0.53     0.19     8.5  -1.5e+01  -0.37   1.3e+01     2008      992     1.00
theta[53]         9.1      1.4      10  -8.3e+00    9.4   2.5e+01       54       27      1.1
theta[54]         -14      1.9      12  -3.2e+01    -13   5.2e+00       38       19      1.1
theta[55]         6.9      1.3      10  -1.0e+01    7.2   2.2e+01       62       31      1.1
theta[56]        -9.0      1.5      11  -2.5e+01   -9.3   8.7e+00       48       24      1.1
theta[57]         6.3     0.62     9.2  -9.0e+00    6.7   2.1e+01      220      109      1.0
theta[58]         9.4     0.69     9.1  -5.7e+00     10   2.4e+01      173       85      1.0
theta[59]         8.0     0.32     8.8  -6.5e+00    8.4   2.3e+01      757      374      1.0
theta[60]          22      2.9      15  -5.6e-01     22   4.5e+01       26       13      1.2
theta[61]        -9.8      2.0      12  -2.9e+01   -9.3   1.0e+01       36       18      1.1
theta[62]         6.1     0.59     8.7  -8.8e+00    6.7   2.0e+01      223      110      1.0
theta[63]         8.0     0.79     8.9  -6.9e+00    8.6   2.2e+01      127       63      1.0
theta[64]        0.10     0.25     8.5  -1.3e+01  -0.46   1.5e+01     1181      583      1.0
theta[65]         7.5     0.77     9.2  -7.6e+00    8.2   2.2e+01      143       71      1.0
Phi              0.36    0.016    0.24  -2.1e-02   0.36   7.6e-01      241      119      1.0
R[1,1]            178       26     120   6.7e+00    177   3.8e+02       22       11      1.2
Q[1,1]            203       24     122   2.8e+01    191   4.2e+02       26       13      1.2
y_hat[1]           56     0.16     9.3   4.0e+01     56   7.1e+01     3285     1622      1.0
y_hat[2]           50     0.55     8.8   3.6e+01     49   6.5e+01      259      128      1.0
y_hat[3]           58     0.22     8.1   4.4e+01     58   7.1e+01     1368      676      1.0
y_hat[4]           59     0.19     8.2   4.6e+01     60   7.3e+01     1792      885      1.0
y_hat[5]           62     0.32     8.6   4.7e+01     62   7.6e+01      724      358      1.0
y_hat[6]           66      1.1     9.7   5.0e+01     67   8.0e+01       78       38      1.1
y_hat[7]           60     0.24     8.2   4.6e+01     60   7.4e+01     1167      576      1.0
y_hat[8]           60     0.98     9.1   4.4e+01     61   7.3e+01       86       42      1.0
y_hat[9]           40      1.5      11   2.4e+01     39   5.7e+01       49       24      1.1
y_hat[10]          40     0.85     9.9   2.6e+01     39   5.7e+01      137       67      1.0
y_hat[11]          36      1.5      11   2.1e+01     35   5.5e+01       54       27      1.1
y_hat[12]          40     0.91     9.7   2.5e+01     39   5.7e+01      114       56      1.0
y_hat[13]          42      1.1     9.7   2.7e+01     41   5.9e+01       84       41      1.1
y_hat[14]          54     0.14     8.0   4.0e+01     54   6.7e+01     3084     1523      1.0
y_hat[15]          63      1.5      11   4.5e+01     64   7.9e+01       54       26      1.1
y_hat[16]          45     0.76     8.9   3.1e+01     44   5.9e+01      138       68      1.0
y_hat[17]          34      2.0      12   1.6e+01     33   5.5e+01       38       19      1.1
y_hat[18]          36      1.9      12   1.9e+01     35   5.6e+01       40       20      1.1
y_hat[19]          55     0.46     8.5   4.0e+01     56   6.8e+01      339      167      1.0
y_hat[20]          52     0.53     8.4   3.9e+01     51   6.6e+01      247      122      1.0
y_hat[21]          65      1.0     9.3   4.9e+01     66   7.8e+01       79       39      1.1
y_hat[22]          65     0.95     9.3   4.9e+01     66   7.9e+01       96       47      1.1
y_hat[23]          55     0.29     8.4   4.3e+01     55   7.0e+01      842      416      1.0
y_hat[24]          64     0.98     9.3   4.9e+01     65   7.8e+01       91       45      1.1
y_hat[25]          59     0.53     8.7   4.4e+01     60   7.2e+01      267      132      1.0
y_hat[26]          44      1.2     9.7   3.0e+01     43   6.0e+01       69       34      1.1
y_hat[27]          47     0.28     8.2   3.3e+01     46   6.1e+01      838      414      1.0
y_hat[28]          39      1.5      11   2.4e+01     38   5.8e+01       54       27      1.1
y_hat[29]          47     0.38     8.4   3.3e+01     47   6.1e+01      484      239      1.0
y_hat[30]          55     0.28     8.2   4.1e+01     56   6.9e+01      853      421      1.0
y_hat[31]          52     0.32     8.4   3.9e+01     51   6.6e+01      667      329      1.0
y_hat[32]          63      1.2     9.4   4.8e+01     64   7.7e+01       65       32      1.1
y_hat[33]          54     0.22     8.0   4.2e+01     54   6.9e+01     1393      688      1.0
y_hat[34]          58     0.13     8.3   4.5e+01     58   7.2e+01     4061     2006      1.0
y_hat[35]          76      2.9      15   5.2e+01     75   9.8e+01       28       14      1.2
y_hat[36]          47     0.96     9.1   3.3e+01     46   6.2e+01       90       44      1.0
y_hat[37]          39      1.2      10   2.5e+01     38   5.7e+01       70       34      1.1
y_hat[38]          32      2.2      13   1.3e+01     32   5.4e+01       35       17      1.1
y_hat[39]          47     0.23     8.3   3.4e+01     47   6.1e+01     1244      614      1.0
y_hat[40]          56     0.69     8.7   4.1e+01     57   6.9e+01      158       78      1.0
y_hat[41]          42      1.5      10   2.7e+01     41   5.9e+01       48       24      1.1
y_hat[42]          48     0.47     8.6   3.4e+01     47   6.2e+01      328      162      1.0
y_hat[43]          53     0.13     8.1   4.0e+01     53   6.7e+01     3742     1848      1.0
y_hat[44]          66     0.80     9.1   5.0e+01     66   8.0e+01      129       64      1.0
y_hat[45]          74      1.6      12   5.4e+01     75   9.0e+01       51       25      1.1
y_hat[46]          70     0.76     9.6   5.3e+01     71   8.4e+01      157       77      1.0
y_hat[47]          72      1.1      10   5.4e+01     73   8.6e+01       86       43      1.1
y_hat[48]          69     0.82     9.7   5.2e+01     71   8.4e+01      140       69      1.0
y_hat[49]          71      1.4      10   5.3e+01     72   8.6e+01       58       29      1.1
y_hat[50]          63     0.66     8.8   4.8e+01     64   7.6e+01      179       88      1.0
y_hat[51]          44      1.8      11   2.8e+01     44   6.2e+01       40       20      1.1
y_hat[52]          55     0.12     8.0   4.1e+01     55   6.8e+01     4664     2303     1.00
y_hat[53]          64      1.4      10   4.7e+01     65   7.9e+01       54       26      1.1
y_hat[54]          41      2.0      12   2.4e+01     42   6.0e+01       36       18      1.1
y_hat[55]          62      1.2     9.8   4.6e+01     62   7.6e+01       62       31      1.1
y_hat[56]          46      1.6      10   3.1e+01     46   6.3e+01       41       20      1.1
y_hat[57]          61     0.55     8.9   4.6e+01     62   7.5e+01      262      129      1.0
y_hat[58]          64     0.66     8.9   4.9e+01     65   7.9e+01      178       88      1.0
y_hat[59]          63     0.28     8.6   4.9e+01     63   7.7e+01      949      468      1.0
y_hat[60]          78      2.9      14   5.5e+01     77   9.9e+01       25       12      1.2
y_hat[61]          45      2.1      12   2.7e+01     46   6.5e+01       34       17      1.1
y_hat[62]          61     0.50     8.5   4.7e+01     62   7.4e+01      291      144      1.0
y_hat[63]          63     0.72     8.6   4.8e+01     64   7.6e+01      141       70      1.0
y_hat[64]          55     0.23     8.2   4.2e+01     54   6.9e+01     1219      602      1.0
y_hat[65]          63     0.74     9.1   4.7e+01     63   7.6e+01      150       74      1.0

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at 
convergence, R_hat=1).

We can find there are sufficient ESS for each parameters and the \(\hat{R}\)’s are close to the 1.00, showing the MCMC chains could be convergent.

y_hat_pos_summary <- ssm0_fit$summary("y_hat", mean, median, quantile2)


s20_pos_predict <- subdata %>% 
  filter(Subject == 20, Affect == "Pos") %>% 
  drop_na(Score) %>% 
  bind_cols(y_hat_pos_summary)

ggplot(s20_pos_predict, aes(x = DateTime, y = Score)) + 
  geom_line(color = pos_neg_color[2]) + 
  geom_point(color = pos_neg_color[2]) +
  geom_line(aes(y = mean), linetype = "dashed") +
  geom_ribbon(aes(ymin = q5, ymax = q95), alpha = 0.25) +
  scale_y_continuous(limits = c(0, 100)) +
  scale_x_datetime(breaks = as_datetime(1:7 * 86400),
                   labels = paste("Day", 1:7),
                   limits = as_datetime(c(1, 8) * 86400))
Figure 2: Postive affect observervation vs. fitted value (the black dotted line) and 95% CI (gray areas).

The fitted line (dotted black line) with 95% CI (gray areas) indeed capture the trend of the observed data.

3.1.4.2 Case 2: Fit PA and NA simultaneously

Here, we fit the PA and NA at the same time

ssm1.stan
#include ssm-function.stan

data {
  int<lower=1> T;
  array[T] vector[2] y;
  vector[2] m_0;
  cov_matrix[2] C_0;
}

parameters {
  vector[2] mu;
  vector[2] theta_0;
  array[T] vector[2] theta;
  matrix[2, 2] Phi;
  cov_matrix[2] R;
  cov_matrix[2] Q;
}

model {
  theta_0 ~ multi_normal(m_0, C_0);
  
  theta[1] ~ multi_normal(Phi * theta_0, Q);
  y[1] ~ multi_normal(mu + theta[1], R);
  
  for (t in 2:T) {
    theta[t] ~ multi_normal(Phi * theta[t-1], Q);
    y[t] ~ multi_normal(mu + theta[t], R);
  }
  
}

generated quantities {
  array[T] vector[2] y_hat;
  matrix[2, 2] Tau; 
  real rel_1; 
  real rel_2;
  
  for (t in 1:T) {
    y_hat[t] = mu + theta[t];
  }
  
  Tau = to_matrix((identity_matrix(4) - kronecker_prod(Phi, Phi)) \ to_vector(Q), 2, 2);
  
  rel_1 = Tau[1, 1] / (Tau[1, 1] + R[1, 1]);
  rel_2 = Tau[2, 2] / (Tau[2, 2] + R[2, 2]);
}
s20_data <- data %>% 
  filter(Subject == 20) %>% 
  drop_na(Score) %>% 
  pivot_wider(names_from = Affect, values_from = Score) %>% 
  select(Pos, Neg)

ssm1 <- cmdstan_model("stan/ssm1.stan")
  
ssm1_data <- list(`T` = 65,
                  y = s20_data, 
                  m_0 = c(0, 0), 
                  C_0 = diag(c(10^3, 10^3), 2, 2))

ssm1_fit <- ssm1$sample(data = ssm1_data, 
                        seed = 20240207,
                        chains = 4,
                        parallel_chains = 4,
                        refresh = 500)

ssm1_fit$save_object(file = "stan/ssm1-fit.RDS")
ssm1_fit <- readRDS("stan/ssm1-fit.RDS")
ssm1_fit$cmdstan_summary()
Inference for Stan model: ssm1_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.

Warmup took (23, 25, 31, 19) seconds, 1.6 minutes total
Sampling took (18, 14, 16, 18) seconds, 1.1 minutes total

                    Mean     MCSE   StdDev     5%       50%    95%  N_Eff  N_Eff/s    R_hat

lp__            -6.6e+02      2.8       30   -706  -6.7e+02   -607    111      1.7      1.0
accept_stat__       0.76  2.3e-02  3.1e-01  0.029      0.92   1.00    176      2.7  1.0e+00
stepsize__         0.054  9.8e-05  4.4e-03  0.047     0.054  0.060   2000       30  1.9e+13
treedepth__          6.1  7.9e-02  9.9e-01    4.0       6.0    7.0    155      2.4  1.0e+00
n_leapfrog__         100  3.2e+00  6.3e+01     30        63    255    401      6.1  1.0e+00
divergent__        0.072  1.9e-02  2.6e-01   0.00      0.00    1.0    190      2.9  1.0e+00
energy__             734  2.9e+00  3.1e+01    679       737    781    115      1.7  1.0e+00

mu[1]            5.5e+01     0.52      5.6     46   5.5e+01     64    120      1.8      1.1
mu[2]            2.5e+01     0.52      5.6     17   2.6e+01     35    115      1.8      1.1
theta_0[1]      -3.5e+00     0.55       24    -44  -3.8e+00     38   1942       30      1.0
theta_0[2]      -2.3e-01     0.74       21    -34  -2.0e-01     33    768       12      1.0
theta[1,1]       1.7e+00     0.62       12    -17   1.4e+00     22    371      5.7      1.0
theta[1,2]       2.4e+00     0.79       10    -15   2.6e+00     19    166      2.5      1.0
theta[2,1]      -4.6e+00     0.66       11    -23  -4.3e+00     13    254      3.9      1.0
theta[2,2]      -3.1e+00     0.56      9.0    -17  -3.6e+00     13    256      3.9      1.0
theta[3,1]       2.9e+00     0.60       10    -14   3.4e+00     20    283      4.3      1.0
theta[3,2]       4.0e+00     0.70      9.4    -12   4.3e+00     19    181      2.8      1.0
theta[4,1]       4.7e-01     0.58      9.6    -16   7.6e-01     15    272      4.1      1.0
theta[4,2]      -2.6e+00     0.59      8.8    -16  -2.7e+00     13    223      3.4      1.0
theta[5,1]       3.8e+00     0.69     10.0    -12   3.5e+00     20    211      3.2      1.0
theta[5,2]       4.8e-01     0.77      8.9    -16   5.5e-01     15    136      2.1      1.0
theta[6,1]       5.1e+00     0.70       11    -12   4.8e+00     23    238      3.6      1.0
theta[6,2]       3.3e+00     0.79      9.5    -14   3.7e+00     19    144      2.2      1.0
theta[7,1]      -5.4e-01     0.75      9.8    -16  -1.2e+00     17    171      2.6      1.0
theta[7,2]       4.1e+00     0.74      8.9    -11   4.5e+00     18    144      2.2      1.0
theta[8,1]      -2.9e-01     0.73       10    -16  -5.4e-01     16    191      2.9      1.0
theta[8,2]       6.2e+00     0.84      9.2  -10.0   6.2e+00     21    120      1.8      1.0
theta[9,1]      -1.4e+01     0.72       11    -33  -1.4e+01    2.2    231      3.5      1.0
theta[9,2]       1.5e+01     0.74      9.7   -1.4   1.5e+01     31    173      2.6      1.0
theta[10,1]     -1.4e+01     0.85       11    -34  -1.4e+01    3.7    177      2.7      1.0
theta[10,2]      1.7e+01     0.88      9.9    1.1   1.7e+01     33    126      1.9      1.0
theta[11,1]     -1.8e+01     0.71       11    -37  -1.7e+01  -0.24    257      3.9      1.0
theta[11,2]      1.8e+01     0.72      9.8    1.5   1.8e+01     34    186      2.8      1.0
theta[12,1]     -1.4e+01     0.68       11    -33  -1.3e+01    2.0    257      3.9      1.0
theta[12,2]      6.8e+00     0.58      8.8   -7.0   6.6e+00     22    230      3.5      1.0
theta[13,1]     -1.0e+01     0.71       11    -28  -9.4e+00    6.1    219      3.3      1.0
theta[13,2]      4.9e+00     0.60      8.9   -8.7   4.7e+00     20    217      3.3      1.0
theta[14,1]     -2.5e+00     0.57      9.7    -20  -2.2e+00     12    285      4.3      1.0
theta[14,2]      4.9e+00     0.56      8.7   -9.0   4.7e+00     19    236      3.6      1.0
theta[15,1]      2.1e+00     0.77       11    -15   2.1e+00     20    204      3.1      1.0
theta[15,2]      3.7e+00     0.65      9.2    -11   3.6e+00     19    201      3.1      1.0
theta[16,1]     -1.2e+01     0.64      9.8    -28  -1.1e+01    2.6    234      3.6      1.0
theta[16,2]      1.4e+01     0.67      8.6   0.20   1.4e+01     29    162      2.5      1.0
theta[17,1]     -2.0e+01     0.93       13    -41  -1.9e+01  -0.62    185      2.8      1.0
theta[17,2]      2.3e+01     0.81       11    6.1   2.3e+01     41    181      2.8      1.0
theta[18,1]     -2.0e+01     0.81       12    -41  -1.9e+01   -1.1    227      3.5      1.0
theta[18,2]      2.4e+01     0.77       11    6.4   2.3e+01     42    202      3.1      1.0
theta[19,1]     -6.4e+00     0.87       11    -24  -6.4e+00     11    156      2.4      1.0
theta[19,2]      1.0e+01     0.58      8.8   -3.4   1.0e+01     25    226      3.4      1.0
theta[20,1]     -6.3e+00     0.87       10    -23  -6.3e+00     12    135      2.1      1.0
theta[20,2]      4.0e+00     0.54      8.4   -9.5   3.8e+00     18    241      3.7      1.0
theta[21,1]      5.1e+00     0.75       11    -12   4.7e+00     24    206      3.1      1.0
theta[21,2]     -3.3e+00     0.60      9.0    -17  -3.6e+00     11    225      3.4      1.0
theta[22,1]      5.6e+00     0.73       11    -11   5.0e+00     23    213      3.2      1.0
theta[22,2]     -3.7e+00     0.63      9.1    -19  -3.9e+00     12    213      3.3      1.0
theta[23,1]     -7.5e-01     0.66       10    -19  -7.3e-01     17    235      3.6      1.0
theta[23,2]      4.9e+00     0.64      9.1    -11   5.2e+00     19    201      3.1      1.0
theta[24,1]      4.2e+00     0.61       10    -12   3.5e+00     22    278      4.2      1.0
theta[24,2]     -1.1e+00     0.59      8.8    -15  -1.5e+00     14    219      3.3      1.0
theta[25,1]      1.4e+00     0.64       10    -16   1.6e+00     19    265      4.0      1.0
theta[25,2]     -5.0e+00     0.63      9.3    -20  -5.2e+00     10    219      3.3      1.0
theta[26,1]     -7.5e+00     0.68       11    -26  -6.8e+00    9.1    261      4.0      1.0
theta[26,2]      9.7e+00     0.70      9.3   -5.5   9.6e+00     25    181      2.8      1.0
theta[27,1]     -8.2e+00     0.52      9.6    -24  -7.7e+00    7.4    342      5.2      1.0
theta[27,2]      8.8e+00     0.50      8.3   -4.3   8.5e+00     23    269      4.1      1.0
theta[28,1]     -1.2e+01     0.69       11    -31  -1.1e+01    4.4    243      3.7      1.0
theta[28,2]      7.8e+00     0.63      8.9   -6.6   7.8e+00     23    200      3.0      1.0
theta[29,1]     -6.1e+00     0.67      9.8    -22  -5.9e+00     11    214      3.3      1.0
theta[29,2]      5.2e-01     0.61      8.6    -13   3.0e-01     15    200      3.1      1.0
theta[30,1]      1.6e+00     0.55      9.5    -14   2.0e+00     17    304      4.6      1.0
theta[30,2]     -3.7e-01     0.55      8.3    -14  -4.1e-01     13    224      3.4      1.0
theta[31,1]     -2.9e+00     0.74       11    -22  -2.5e+00     15    216      3.3      1.0
theta[31,2]     -6.9e+00     0.65      9.4    -22  -7.2e+00    8.8    206      3.1      1.0
theta[32,1]      7.3e+00     0.85       12    -14   7.3e+00     28    213      3.2      1.0
theta[32,2]      6.6e+00     0.79       11    -11   6.7e+00     23    187      2.8      1.0
theta[33,1]     -2.7e+00     0.93      9.8    -18  -3.1e+00     15    109      1.7      1.0
theta[33,2]     -2.2e-01     0.62      8.2    -13  -3.3e-01     14    173      2.6      1.0
theta[34,1]      1.4e+00     0.60      9.1    -13   1.3e+00     17    229      3.5      1.0
theta[34,2]     -2.8e+00     0.54      8.3    -16  -3.0e+00     11    230      3.5      1.0
theta[35,1]      1.3e+01     0.93       14   -7.0   1.2e+01     37    215      3.3      1.0
theta[35,2]     -1.7e+00     0.78       10    -18  -1.5e+00     15    181      2.8      1.0
theta[36,1]     -7.7e+00     0.68      9.9    -24  -7.1e+00    7.7    213      3.2      1.0
theta[36,2]      6.2e+00     0.61      8.6   -7.6   6.5e+00     20    199      3.0      1.0
theta[37,1]     -1.3e+01     0.70       11    -33  -1.3e+01    3.1    246      3.7      1.0
theta[37,2]      8.9e+00     0.58      9.0   -5.5   8.8e+00     24    243      3.7      1.0
theta[38,1]     -1.7e+01      1.2       16    -45  -1.6e+01    6.3    180      2.7      1.0
theta[38,2]      2.4e+01      1.1       14   0.93   2.4e+01     46    149      2.3      1.0
theta[39,1]     -1.0e+01     0.70       10    -28  -1.0e+01    6.8    222      3.4      1.0
theta[39,2]      6.7e+00     0.61      8.7   -7.2   6.7e+00     22    201      3.1      1.0
theta[40,1]     -6.7e-01     0.61       10    -17  -6.5e-01     16    277      4.2      1.0
theta[40,2]      1.5e+00     0.59      8.8    -14   1.4e+00     16    219      3.3      1.0
theta[41,1]     -8.5e+00     0.76       12    -28  -7.8e+00    8.9    229      3.5      1.0
theta[41,2]      1.3e+00     0.63      9.2    -14   1.0e+00     17    217      3.3      1.0
theta[42,1]     -1.7e+00     0.73       11    -20  -1.4e+00     16    224      3.4      1.0
theta[42,2]     -6.0e+00     0.64      9.3    -21  -6.1e+00    9.3    212      3.2      1.0
theta[43,1]      3.0e+00     0.62       10    -15   3.2e+00     19    277      4.2      1.0
theta[43,2]     -1.2e+01     0.62      9.1    -26  -1.2e+01    3.3    218      3.3      1.0
theta[44,1]      1.1e+01     0.57      9.8   -4.1   1.1e+01     27    297      4.5      1.0
theta[44,2]     -1.1e+01     0.51      8.4    -25  -1.1e+01    2.6    271      4.1      1.0
theta[45,1]      1.6e+01     0.70       11   -1.4   1.5e+01     34    252      3.8      1.0
theta[45,2]     -1.0e+01     0.61      9.1    -25  -1.0e+01    4.6    225      3.4      1.0
theta[46,1]      1.1e+01     0.69       10   -5.2   1.1e+01     29    225      3.4      1.0
theta[46,2]     -8.2e+00     0.60      8.9    -24  -8.2e+00    6.3    221      3.4      1.0
theta[47,1]      1.3e+01     0.66       10   -3.4   1.2e+01     30    250      3.8      1.0
theta[47,2]     -8.3e+00     0.53      8.7    -23  -8.2e+00    6.0    275      4.2      1.0
theta[48,1]      1.2e+01     0.81       11   -5.5   1.1e+01     32    183      2.8      1.0
theta[48,2]     -1.1e+01     0.63      8.9    -25  -1.1e+01    3.7    199      3.0      1.0
theta[49,1]      1.4e+01     0.77       11   -3.6   1.3e+01     33    198      3.0      1.0
theta[49,2]     -1.2e+01     0.65      9.1    -27  -1.2e+01    3.0    196      3.0      1.0
theta[50,1]      7.9e+00     0.57      9.4   -6.9   7.6e+00     24    269      4.1      1.0
theta[50,2]     -7.9e+00     0.56      8.5    -22  -7.9e+00    6.5    225      3.4      1.0
theta[51,1]     -4.1e+00     0.79       12    -24  -3.6e+00     14    213      3.2      1.0
theta[51,2]     -2.1e+00     0.68      9.6    -18  -2.3e+00     14    200      3.0      1.0
theta[52,1]      2.5e+00     0.54      9.7    -14   2.7e+00     18    318      4.8      1.0
theta[52,2]     -1.0e+01     0.59      8.8    -24  -1.0e+01    4.7    221      3.4      1.0
theta[53,1]      1.1e+01     0.60      9.9   -5.8   1.1e+01     27    274      4.2      1.0
theta[53,2]     -8.9e+00     0.69      8.7    -22  -9.0e+00    5.7    159      2.4      1.0
theta[54,1]     -5.7e+00     0.78       12    -28  -4.4e+00     13    251      3.8      1.0
theta[54,2]     -4.7e+00     0.69     10.0    -20  -5.3e+00     13    209      3.2      1.0
theta[55,1]      9.9e+00     0.56       10   -6.7   1.0e+01     26    335      5.1      1.0
theta[55,2]     -9.8e+00     0.57      8.7    -23  -1.0e+01    4.3    230      3.5      1.0
theta[56,1]     -2.2e+00     0.84       13    -24  -1.3e+00     17    223      3.4      1.0
theta[56,2]     -1.2e+01     0.70       10    -29  -1.2e+01    4.9    220      3.4      1.0
theta[57,1]      9.5e+00     0.48      9.9   -6.2   9.3e+00     26    413      6.3      1.0
theta[57,2]     -7.7e+00     0.45      8.3    -21  -7.6e+00    5.8    346      5.3      1.0
theta[58,1]      9.1e+00     0.48      9.5   -5.9   8.3e+00     26    400      6.1      1.0
theta[58,2]     -7.7e+00     0.46      8.1    -21  -7.9e+00    5.4    306      4.7      1.0
theta[59,1]      7.3e+00     0.49      9.4   -7.7   7.0e+00     24    367      5.6      1.0
theta[59,2]     -1.1e+01     0.50      8.3    -25  -1.1e+01    2.4    274      4.2      1.0
theta[60,1]      1.8e+01     0.77       12  -0.31   1.7e+01     39    253      3.9      1.0
theta[60,2]     -1.0e+01     0.63      9.6    -26  -1.0e+01    5.7    234      3.6      1.0
theta[61,1]     -4.9e+00     0.71       12    -26  -3.8e+00     12    266      4.1      1.0
theta[61,2]     -5.6e-01     0.56      9.3    -16  -7.2e-01     15    270      4.1      1.0
theta[62,1]      5.6e+00     0.51      9.3    -10   5.6e+00     21    340      5.2      1.0
theta[62,2]     -6.3e+00     0.51      8.3    -20  -6.4e+00    7.4    267      4.1      1.0
theta[63,1]      5.7e+00     0.60      9.9    -10   5.3e+00     23    272      4.1      1.0
theta[63,2]     -7.3e+00     0.66      8.5    -21  -7.5e+00    7.3    166      2.5      1.0
theta[64,1]      1.4e+00     0.74       12    -19   1.7e+00     21    261      4.0      1.0
theta[64,2]      9.1e+00     0.87       11   -9.3   8.8e+00     28    162      2.5      1.0
theta[65,1]      3.2e+00     0.83       12    -15   3.2e+00     22    194      3.0      1.0
theta[65,2]     -7.4e+00     0.68      9.4    -23  -7.6e+00    7.8    192      2.9      1.0
Phi[1,1]        -6.0e-03    0.037     0.52  -0.84  -2.0e-03   0.87    197      3.0      1.0
Phi[1,2]        -6.4e-01    0.037     0.56   -1.6  -6.0e-01   0.19    229      3.5      1.0
Phi[2,1]        -6.4e-02    0.038     0.54  -0.93  -7.3e-02   0.77    205      3.1      1.0
Phi[2,2]         5.8e-01    0.036     0.52  -0.26   5.7e-01    1.4    204      3.1      1.0
R[1,1]           2.5e+02      6.8      103     78   2.5e+02    420    233      3.6      1.0
R[1,2]          -1.2e+02      4.8       69   -236  -1.2e+02    -14    209      3.2      1.0
R[2,1]          -1.2e+02      4.8       69   -236  -1.2e+02    -14    209      3.2      1.0
R[2,2]           1.3e+02      4.3       59     33   1.2e+02    224    190      2.9      1.0
Q[1,1]           1.3e+02      7.1       95     15   1.0e+02    311    178      2.7      1.0
Q[1,2]          -7.4e+01      5.2       68   -205  -6.0e+01     12    168      2.6      1.0
Q[2,1]          -7.4e+01      5.2       68   -205  -6.0e+01     12    168      2.6      1.0
Q[2,2]           1.0e+02      5.2       65     18   9.3e+01    225    152      2.3      1.0
y_hat[1,1]       5.7e+01     0.31       11     40   5.7e+01     75   1188       18      1.0
y_hat[1,2]       2.8e+01     0.29      8.1     15   2.8e+01     40    762       12      1.0
y_hat[2,1]       5.1e+01     0.38      9.4     35   5.1e+01     67    594      9.1      1.0
y_hat[2,2]       2.2e+01     0.24      7.3     10   2.2e+01     35    929       14      1.0
y_hat[3,1]       5.8e+01     0.34      8.8     44   5.8e+01     73    689       11      1.0
y_hat[3,2]       2.9e+01     0.35      7.3     17   3.0e+01     41    435      6.6      1.0
y_hat[4,1]       5.6e+01     0.29      8.3     42   5.6e+01     70    808       12      1.0
y_hat[4,2]       2.3e+01     0.29      7.1     11   2.3e+01     35    595      9.1      1.0
y_hat[5,1]       5.9e+01     0.28      8.5     46   5.9e+01     74    952       15      1.0
y_hat[5,2]       2.6e+01     0.27      7.0     14   2.6e+01     37    685       10      1.0
y_hat[6,1]       6.1e+01     0.46      9.8     45   6.0e+01     77    466      7.1      1.0
y_hat[6,2]       2.9e+01     0.39      7.8     15   2.9e+01     41    408      6.2      1.0
y_hat[7,1]       5.5e+01     0.31      8.6     41   5.5e+01     69    746       11      1.0
y_hat[7,2]       3.0e+01     0.30      7.0     18   3.0e+01     41    557      8.5      1.0
y_hat[8,1]       5.5e+01     0.36      8.8     41   5.5e+01     70    594      9.1      1.0
y_hat[8,2]       3.2e+01     0.31      7.3     19   3.2e+01     43    568      8.7      1.0
y_hat[9,1]       4.1e+01     0.46      9.7     25   4.1e+01     56    441      6.7      1.0
y_hat[9,2]       4.0e+01     0.42      7.8     27   4.1e+01     53    349      5.3      1.0
y_hat[10,1]      4.1e+01     0.50      9.6     25   4.1e+01     57    377      5.7      1.0
y_hat[10,2]      4.3e+01     0.52      7.9     30   4.3e+01     55    228      3.5      1.0
y_hat[11,1]      3.8e+01     0.52       10     21   3.8e+01     54    382      5.8      1.0
y_hat[11,2]      4.3e+01     0.44      8.0     31   4.4e+01     57    334      5.1      1.0
y_hat[12,1]      4.1e+01     0.41      9.6     25   4.2e+01     57    546      8.3      1.0
y_hat[12,2]      3.2e+01     0.30      7.3     20   3.2e+01     44    588      9.0      1.0
y_hat[13,1]      4.5e+01     0.38      9.2     30   4.6e+01     60    569      8.7      1.0
y_hat[13,2]      3.0e+01     0.22      7.0     19   3.0e+01     42   1001       15      1.0
y_hat[14,1]      5.3e+01     0.30      8.4     39   5.3e+01     66    802       12      1.0
y_hat[14,2]      3.0e+01     0.21      6.7     19   3.1e+01     41   1019       16      1.0
y_hat[15,1]      5.8e+01     0.53      9.9     42   5.7e+01     75    346      5.3      1.0
y_hat[15,2]      2.9e+01     0.36      7.5     16   2.9e+01     41    441      6.7      1.0
y_hat[16,1]      4.4e+01     0.33      8.5     30   4.4e+01     58    670       10      1.0
y_hat[16,2]      4.0e+01     0.19      6.7     28   4.0e+01     51   1207       18      1.0
y_hat[17,1]      3.6e+01     0.73       11     16   3.6e+01     53    248      3.8      1.0
y_hat[17,2]      4.9e+01     0.59      9.3     34   4.8e+01     63    244      3.7      1.0
y_hat[18,1]      3.5e+01     0.65       11     16   3.6e+01     54    301      4.6      1.0
y_hat[18,2]      4.9e+01     0.58      9.5     34   4.9e+01     65    266      4.1      1.0
y_hat[19,1]      4.9e+01     0.46      9.4     34   4.9e+01     64    419      6.4      1.0
y_hat[19,2]      3.6e+01     0.24      6.9     24   3.6e+01     47    796       12      1.0
y_hat[20,1]      4.9e+01     0.38      8.4     35   4.9e+01     63    499      7.6      1.0
y_hat[20,2]      2.9e+01     0.21      6.5     19   3.0e+01     39    987       15      1.0
y_hat[21,1]      6.1e+01     0.48      9.3     46   6.0e+01     75    378      5.8      1.0
y_hat[21,2]      2.2e+01     0.35      7.3    9.9   2.2e+01     34    436      6.7      1.0
y_hat[22,1]      6.1e+01     0.54      9.7     47   6.1e+01     77    320      4.9      1.0
y_hat[22,2]      2.2e+01     0.39      7.5    9.4   2.2e+01     34    374      5.7      1.0
y_hat[23,1]      5.5e+01     0.31      8.6     40   5.5e+01     68    771       12      1.0
y_hat[23,2]      3.0e+01     0.38      7.2     18   3.0e+01     42    370      5.6      1.0
y_hat[24,1]      6.0e+01     0.46      9.3     45   5.9e+01     76    412      6.3      1.0
y_hat[24,2]      2.4e+01     0.32      7.3     13   2.4e+01     36    516      7.9      1.0
y_hat[25,1]      5.7e+01     0.50      9.6     42   5.7e+01     73    374      5.7      1.0
y_hat[25,2]      2.0e+01     0.43      8.1    7.0   2.0e+01     34    349      5.3      1.0
y_hat[26,1]      4.8e+01     0.50      9.7     32   4.8e+01     63    381      5.8      1.0
y_hat[26,2]      3.5e+01     0.43      7.5     23   3.5e+01     47    304      4.6      1.0
y_hat[27,1]      4.7e+01     0.43      8.9     32   4.8e+01     62    424      6.5      1.0
y_hat[27,2]      3.4e+01     0.24      6.8     23   3.4e+01     45    765       12      1.0
y_hat[28,1]      4.4e+01     0.50      9.5     28   4.4e+01     59    362      5.5      1.0
y_hat[28,2]      3.3e+01     0.36      7.3     21   3.3e+01     45    425      6.5      1.0
y_hat[29,1]      4.9e+01     0.38      8.7     35   5.0e+01     63    508      7.7      1.0
y_hat[29,2]      2.6e+01     0.33      7.0     15   2.5e+01     38    452      6.9      1.0
y_hat[30,1]      5.7e+01     0.23      8.0     44   5.7e+01     70   1258       19      1.0
y_hat[30,2]      2.5e+01     0.15      6.3     15   2.5e+01     35   1748       27      1.0
y_hat[31,1]      5.3e+01     0.45      9.6     36   5.3e+01     68    448      6.8      1.0
y_hat[31,2]      1.8e+01     0.42      7.9    5.6   1.9e+01     31    359      5.5      1.0
y_hat[32,1]      6.3e+01     0.65       11     45   6.2e+01     82    298      4.5      1.0
y_hat[32,2]      3.2e+01     0.53      9.0     17   3.2e+01     46    290      4.4      1.0
y_hat[33,1]      5.3e+01     0.30      8.2     39   5.3e+01     65    753       11      1.0
y_hat[33,2]      2.5e+01     0.16      6.2     15   2.5e+01     36   1488       23      1.0
y_hat[34,1]      5.7e+01     0.27      7.8     44   5.7e+01     69    834       13      1.0
y_hat[34,2]      2.3e+01     0.22      6.6     11   2.3e+01     34    928       14      1.0
y_hat[35,1]      6.9e+01     0.77       13     50   6.7e+01     92    271      4.1      1.0
y_hat[35,2]      2.4e+01     0.50      9.0    8.8   2.4e+01     38    322      4.9      1.0
y_hat[36,1]      4.8e+01     0.34      8.6     33   4.8e+01     62    646      9.8      1.0
y_hat[36,2]      3.2e+01     0.25      6.9     20   3.2e+01     43    767       12      1.0
y_hat[37,1]      4.2e+01     0.46      9.8     25   4.2e+01     58    454      6.9      1.0
y_hat[37,2]      3.4e+01     0.33      7.6     22   3.4e+01     46    536      8.2      1.0
y_hat[38,1]      3.8e+01      1.0       15     12   3.9e+01     61    214      3.3      1.0
y_hat[38,2]      4.9e+01     0.95       12     30   4.9e+01     70    168      2.6      1.0
y_hat[39,1]      4.5e+01     0.28      8.8     30   4.5e+01     59    969       15      1.0
y_hat[39,2]      3.2e+01     0.19      6.7     22   3.2e+01     44   1308       20      1.0
y_hat[40,1]      5.5e+01     0.35      8.7     40   5.5e+01     68    614      9.4      1.0
y_hat[40,2]      2.7e+01     0.29      6.9     15   2.7e+01     38    562      8.6      1.0
y_hat[41,1]      4.7e+01     0.51     10.0     30   4.8e+01     62    389      5.9      1.0
y_hat[41,2]      2.7e+01     0.35      7.4     15   2.7e+01     39    464      7.1      1.0
y_hat[42,1]      5.4e+01     0.53      9.5     38   5.4e+01     69    329      5.0      1.0
y_hat[42,2]      1.9e+01     0.39      7.7    7.7   1.9e+01     33    397      6.0      1.0
y_hat[43,1]      5.8e+01     0.51      9.4     43   5.9e+01     74    341      5.2      1.0
y_hat[43,2]      1.4e+01     0.53      8.1    1.4   1.3e+01     29    231      3.5      1.0
y_hat[44,1]      6.7e+01     0.28      8.6     53   6.6e+01     81    957       15     1.00
y_hat[44,2]      1.5e+01     0.25      7.0    3.5   1.4e+01     26    794       12      1.0
y_hat[45,1]      7.1e+01     0.48      9.9     55   7.1e+01     88    425      6.5      1.0
y_hat[45,2]      1.5e+01     0.34      7.4    3.5   1.5e+01     28    478      7.3      1.0
y_hat[46,1]      6.7e+01     0.35      8.9     52   6.7e+01     81    641      9.8      1.0
y_hat[46,2]      1.7e+01     0.29      7.1    4.9   1.7e+01     29    601      9.2      1.0
y_hat[47,1]      6.8e+01     0.38      9.1     54   6.8e+01     83    568      8.7      1.0
y_hat[47,2]      1.7e+01     0.26      7.0    5.5   1.7e+01     29    740       11      1.0
y_hat[48,1]      6.7e+01     0.52      9.5     52   6.7e+01     83    339      5.2      1.0
y_hat[48,2]      1.5e+01     0.34      7.3    2.4   1.5e+01     26    445      6.8      1.0
y_hat[49,1]      6.9e+01     0.52      9.6     54   6.9e+01     85    338      5.2      1.0
y_hat[49,2]      1.4e+01     0.35      7.3    1.6   1.3e+01     25    438      6.7      1.0
y_hat[50,1]      6.3e+01     0.29      8.2     50   6.3e+01     77    813       12      1.0
y_hat[50,2]      1.7e+01     0.22      6.8    6.6   1.7e+01     29    913       14      1.0
y_hat[51,1]      5.1e+01     0.65       11     33   5.2e+01     67    270      4.1      1.0
y_hat[51,2]      2.3e+01     0.47      8.2     11   2.3e+01     37    299      4.6      1.0
y_hat[52,1]      5.8e+01     0.41      8.8     43   5.8e+01     72    473      7.2      1.0
y_hat[52,2]      1.5e+01     0.35      7.4    3.9   1.5e+01     28    449      6.8      1.0
y_hat[53,1]      6.6e+01     0.35      8.7     52   6.6e+01     81    607      9.3      1.0
y_hat[53,2]      1.6e+01     0.31      7.0    4.8   1.6e+01     28    524      8.0      1.0
y_hat[54,1]      5.0e+01     0.69       12     30   5.1e+01     67    282      4.3      1.0
y_hat[54,2]      2.1e+01     0.49      8.7    7.2   2.0e+01     36    315      4.8      1.0
y_hat[55,1]      6.5e+01     0.27      8.8     51   6.5e+01     80   1047       16      1.0
y_hat[55,2]      1.6e+01     0.24      7.1    4.1   1.5e+01     27    890       14      1.0
y_hat[56,1]      5.3e+01     0.77       12     32   5.4e+01     72    239      3.6      1.0
y_hat[56,2]      1.3e+01     0.59      9.5   -1.6   1.3e+01     30    262      4.0      1.0
y_hat[57,1]      6.5e+01     0.48      9.1     49   6.5e+01     79    369      5.6      1.0
y_hat[57,2]      1.8e+01     0.23      6.9    6.3   1.8e+01     29    924       14      1.0
y_hat[58,1]      6.5e+01     0.38      8.9     50   6.5e+01     79    544      8.3      1.0
y_hat[58,2]      1.8e+01     0.24      6.8    6.6   1.8e+01     29    824       13      1.0
y_hat[59,1]      6.3e+01     0.39      8.8     48   6.3e+01     78    509      7.8      1.0
y_hat[59,2]      1.4e+01     0.43      7.4    2.8   1.4e+01     26    298      4.5      1.0
y_hat[60,1]      7.3e+01     0.66       11     57   7.2e+01     94    299      4.6      1.0
y_hat[60,2]      1.5e+01     0.44      8.3    1.5   1.6e+01     28    350      5.3      1.0
y_hat[61,1]      5.1e+01     0.55       11     31   5.1e+01     66    365      5.6      1.0
y_hat[61,2]      2.5e+01     0.39      7.8     12   2.5e+01     38    405      6.2      1.0
y_hat[62,1]      6.1e+01     0.37      8.7     48   6.1e+01     75    557      8.5      1.0
y_hat[62,2]      1.9e+01     0.32      7.2    7.5   1.9e+01     30    501      7.6      1.0
y_hat[63,1]      6.1e+01     0.34      8.9     47   6.1e+01     76    675       10      1.0
y_hat[63,2]      1.8e+01     0.30      7.0    6.6   1.8e+01     30    558      8.5      1.0
y_hat[64,1]      5.7e+01     0.59       11     38   5.7e+01     75    363      5.5      1.0
y_hat[64,2]      3.4e+01     0.68      9.8     19   3.4e+01     51    206      3.1      1.0
y_hat[65,1]      5.9e+01     0.55       10     42   5.9e+01     76    366      5.6      1.0
y_hat[65,2]      1.8e+01     0.39      8.2    4.0   1.8e+01     32    446      6.8      1.0
Tau[1,1]         2.6e+02       13      679     51   2.2e+02    544   2914       44     1.00
Tau[1,2]        -1.9e+02       11      660   -430  -1.5e+02    -22   3649       56      1.0
Tau[2,1]        -1.9e+02       11      660   -430  -1.5e+02    -22   3649       56      1.0
Tau[2,2]         2.2e+02       14      860     63   1.8e+02    463   3803       58      1.0
rel_1            5.0e-01    0.017     0.53   0.14   4.8e-01   0.86   1026       16      1.0
rel_2            6.1e-01    0.015     0.33   0.27   6.1e-01   0.92    501      7.6      1.0

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at 
convergence, R_hat=1).

It seems have enough ESS for each parameters and \(\hat{R}\)’s are closed to 1.

y_hat_summary <- ssm1_fit$summary("y_hat", mean, median, quantile2) %>% 
  mutate(Index = str_extract(variable, "\\d+") %>% as.integer(),
         Affect = str_detect(variable, ",1]") %>% ifelse("Pos", "Neg")) 


s20_predict <- subdata %>% 
  filter(Subject == 20) %>% 
  drop_na(Score) %>% 

  mutate(Index = rep(1:(n()/2), times = 2)) %>% 
  left_join(y_hat_summary)
  

ggplot(s20_predict, aes(x = DateTime, y = Score)) + 
  geom_line(aes(color = Affect)) + 
  geom_point(aes(color = Affect)) +
  scale_color_manual(values = pos_neg_color) +
  geom_line(aes(y = mean, group = Affect), linetype = "dashed") +
  geom_ribbon(aes(ymin = q5, ymax = q95, group = Affect), alpha = 0.25) +
  scale_y_continuous(limits = c(0, 100)) +
  scale_x_datetime(breaks = as_datetime(1:7 * 86400),
                   labels = paste("Day", 1:7),
                   limits = as_datetime(c(1, 8) * 86400))

The fitted lines for both PA and RA capture the observed values for Subject 20.

The posterior distributions of the reliability for positive and negative affect measurement are presented following.

ssm1_draws <- ssm1_fit$draws(format = "df") %>% 
  mutate(rel_diff = rel_1 - rel_2)

mcmc_areas(ssm1_draws, 
           pars = c("rel_1", "rel_2", "rel_diff"), 
           prob = 0.8,
           prob_outer = 0.99)
Figure 3: The posterior distributions of the reliability for PA and NA and the different between of them (\(R^{PA} - R^{NA}\))

3.2 Multilevel extension

3.2.1 Model Specification

Level 2 (between subject)

\[ \boldsymbol{\mu} = \begin{pmatrix}\mu_{1 i} \\\mu_{2 i} \end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_\mu = \begin{pmatrix} \gamma_{\mu 1} \\\gamma_{\mu 2} \end{pmatrix}, \boldsymbol{\Psi}_\mu =\begin{pmatrix} \psi_{\mu1}^2 & \psi_{\mu12} \\ \psi_{\mu12} & \psi_{\mu2}^2\end{pmatrix}\right) \]

\[ vec(\boldsymbol{\Phi}_i) = \begin{pmatrix} \phi_{11i} \\ \phi_{12i} \\ \phi_{21i} \\ \phi_{22i} \end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_{\Phi} = \begin{pmatrix}\gamma_{\phi 11} \\ \gamma_{\phi 12} \\ \gamma_{\phi 21} \\ \gamma_{\phi 22} \end{pmatrix}, \boldsymbol{\boldsymbol{\Psi}}_{\phi} = \begin{pmatrix} \psi_{\phi11}^2 & \psi_{\phi11\phi12} & \psi_{\phi11\phi21} & \psi_{\phi11\phi22} \\ \psi_{\phi11\phi12} & \psi_{\phi12}^2 & \psi_{\phi12\phi21} & \psi_{\phi12\phi22} \\ \psi_{\phi21\phi11} & \psi_{\phi21\phi12} & \psi_{\phi21}^2 & \psi_{\phi21\phi22} \\ \psi_{\phi11\phi22} & \psi_{\phi22\phi12} & \psi_{\phi22\phi21} & \psi_{\phi22}^2 \end{pmatrix} \right) \]

\[ vec(upp.tri(\mathbf{R}_i)) = \begin{pmatrix}\sigma_{\epsilon_{11 i}}^{2} \\\sigma_{\epsilon_{12 i}} \\\sigma_{\epsilon_{22 i}}^{2}\end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_{\sigma_\epsilon^2} = \begin{pmatrix} \gamma_{\sigma_{\epsilon 11}^{2}} \\\gamma_{\sigma_{\epsilon 12}} \\ \gamma_{\sigma_{\epsilon 22}^{2}} \end{pmatrix}, \boldsymbol{\Psi}_{\boldsymbol{\sigma}_{\epsilon}^{2}} = \begin{pmatrix} \psi_{\sigma_{\epsilon 11}^2} & & 0 \\ & \psi_{\sigma_{\epsilon 12}} & \\ 0 & & \psi_{\sigma_{\epsilon 22}^2}\end{pmatrix}\right) \]

\[ vec(upp.tri(\mathbf{Q}_i)) = \begin{pmatrix}\sigma_{\omega_{11 i}}^{2} \\\sigma_{\omega_{12 i}} \\\sigma_{\omega_{22 i}}^{2}\end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_{\sigma_\omega^2} = \begin{pmatrix}\gamma_{\sigma_{\omega 11}^{2}} \\\gamma_{\sigma_{\omega 12}} \\\gamma_{\sigma_{\omega 22}^{2}} \end{pmatrix}, \boldsymbol{\Psi}_{\boldsymbol{\sigma}_{\omega}^{2}} = \begin{pmatrix} \psi_{\sigma_{\omega 11}^2} & & 0 \\ & \psi_{\sigma_{\omega 12}} & \\ 0 & & \psi_{\sigma_{\omega 22}^2}\end{pmatrix}\right) \]

However, the possible range of \(\sigma_{\epsilon 11i}^2\) is \((0, \infty)\) and the one of \(\sigma_{\epsilon 12i}\) is \((-\infty, \infty)\). Thus, the

3.2.2 Reliability

Again, according to Schuurman & Hamaker (2019), the reliability for the multilevel MEVAR(1) model have two types:

  • Reliability for systematic between-subject difference

\[ R_k^B = \frac{(\Psi^2_\mu)_{kk}}{Var(\mathbf{y}_k)} = \frac{(\Psi^2_\mu)_{kk}}{(\Psi^2_\mu)_{kk} + E_{i}[(\mathbf{T}_i)_{kk}]+ (\boldsymbol{\gamma}_{\sigma_\epsilon^2})_k} \]

  • Reliability for within-subject fluctuations

\[ R_{ik}^W = \frac{(\mathbf{T}_i)_{kk}}{Var(y_{ik})} = \frac{(\mathbf{T}_i)_{kk}}{(\mathbf{T}_i)_{kk} + (\mathbf{R}_i)_{kk}} \]

3.2.3 Bayesian Estimation

multilevel-ssm.stan
#include ssm-function.stan

data {
  int<lower=1> N; // number of subjects
  array[N] int<lower=1> T; // number of observation for each subject
  int<lower=1> max_T; // maximum number of observation
  int<lower=1> P; // number of affects
  array[N] matrix[P, max_T] y; // observations 
  array[N] vector[P] m_0; // prior mean of the intial state
  array[N] cov_matrix[P] C_0; // prior covariance of the intial state
}

parameters {
  array[N] vector[P] mu; // ground mean/trane
  array[N] vector[P] theta_0; // initial latent state
  array[N] matrix[P, max_T] theta; // latent states
  array[N] matrix[P, P] Phi; // autoregressive parameters
  array[N] cov_matrix[P] R; // covariance of the measurment error
  array[N] cov_matrix[P] Q; // covariance of the innovation noise
  
  vector[P] gamma_mu; // prior mean of the ground mean
  cov_matrix[P] Psi_mu; // prior covariance of the ground mean
  vector[P * P] gamma_Phi; // prior mean of the autoregressive parameters
  cov_matrix[P * P] Psi_Phi; // prior covariance of the autoregressive parameters
  vector[P * (P + 1) / 2] gamma_R; // prior mean of the covariance of the measurement error
  vector[P * (P + 1) / 2] diag_Psi_R; // prior covariance of the covariance of the measurement error
  vector[P * (P + 1) / 2] gamma_Q; // prior mean of the covariance of the innovation noise
  vector[P * (P + 1) / 2] diag_Psi_Q; // prior covariance of the covariance of innovation noise
}

model {
  // level 1 (within subject)
  for (n in 1:N) {
    // when t = 0
    theta_0[n] ~ multi_normal(m_0[n], C_0[n]);
  
    // when t = 1
    theta[n][, 1] ~ multi_normal(Phi[n] * theta_0[n], Q[n]);
    y[n][, 1] ~ multi_normal(mu[n] + theta[n][, 1], R[n]);
    
    // when t = 2, ..., T 
    for (t in 2:T[n]) {
      theta[n][, t] ~ multi_normal(Phi[n] * theta[n][, t - 1], Q[n]);
      y[n][, t] ~ multi_normal(mu[n] + theta[n][, t], R[n]);
    }
  }
  
  // level 2 (between subject)
  for (n in 1:N) {
    mu[n] ~ multi_normal(gamma_mu, Psi_mu);
    to_vector(Phi[n]) ~ multi_normal(gamma_Phi, Psi_Phi);
    to_vector_lower_tri(R[n]) ~ normal(gamma_R, sqrt(diag_Psi_R));
    to_vector_lower_tri(Q[n]) ~ normal(gamma_Q, sqrt(diag_Psi_Q));
  }
  
  // the (hyper)priors of parameters are set as the Stan default values
}

generated quantities {
  array[N] matrix[P, max_T] y_hat;
  array[N] matrix[P, P] Tau; 
  array[N] vector[P] rel_W;
  vector[P] rel_B;
  
  for (n in 1:N) {
    // prediction 
    for (t in 1:T[n]) {
      y_hat[n][, t] = mu[n] + theta[n][, t];
    }
    
    // within-subject reliability
    Tau[n] = to_matrix((identity_matrix(P * P) - kronecker_prod(Phi[n], Phi[n])) \ to_vector(Q[n]), P, P);
  
    for (p in 1:P) {
      rel_W[n, p] = Tau[n, p, p] / (Tau[n, p, p] + R[n, p, p]);
    }
  }
  
  // between-subject reliability
  for (p in 1:P) {
    rel_B[p] = Psi_mu[p, p] / (Psi_mu[p, p] + mean(Tau[, p, p]) + gamma_R[index_of_diag_lower_tri(p, P)]);
  }
}
data_list <- data %>% 
  group_by(Subject) %>% 
  pivot_wider(names_from = Affect, values_from = Score) %>% 
  select(Pos, Neg) %>% 
  drop_na(Pos, Neg) %>% 
  nest() %>% ungroup() %>% 
  mutate(`T` = map_dbl(data, nrow),
         max_T = max(`T`),
         data_padding = pmap(list(data, `T`, max_T), 
                             \(x, y, z) {
                               bind_rows(x, 
                                         tibble(Pos = rep(Inf, z - y),
                                                Neg = rep(Inf, z - y))) %>% 
                                 t()
                             }))


mssm <- cmdstan_model("stan/multilevel-ssm.stan")

mssm_data <- lst(N = nrow(data_list),
                 `T` = map(data_list$data, nrow),
                 max_T = max(data_list$`T`),
                 P = 2,
                 y = data_list$data_padding, 
                 m_0 = rep(list(c(0, 0)), N), 
                 C_0 = rep(list(diag(c(10^3, 10^3), 2, 2)), N))

mssm_fit <- mssm$sample(data = mssm_data, 
                        seed = 20240221,
                        chains = 4,
                        parallel_chains = 4,
                        iter_sampling = 3000,
                        refresh = 1000, 
                        show_messages = FALSE)
mssm_fit$save_object(file = "stan/multilevel-ssm-fit.RDS")

It took me more than 15 hours to finish the sampling …

mssm_fit <- readRDS("stan/multilevel-ssm-fit.RDS")
mssm_fit_summary <- readRDS("stan/multilevel-ssm-summary.RDS")

skimr::skim(data.frame(Rhat = as.double(mssm_fit_summary$rhat),
                       ESS = as.double(mssm_fit_summary$ess_bulk)))
Data summary
Name data.frame(Rhat = as.doub…
Number of rows 30499
Number of columns 2
_______________________
Column type frequency:
numeric 2
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Rhat 2522 0.92 1.27 0.40 1.01 1.07 1.12 1.24 3.93 ▇▁▁▁▁
ESS 2522 0.92 39.90 45.39 4.30 12.25 23.14 48.24 562.51 ▇▁▁▁▁

However, the most of the ESS are super low. And \(\hat{R}\)’s are deviated form 1, meaning the MCMC samples are not converged yet. Thus, the posterior results are not reliable.

y_hat_summary <- mssm_fit$summary("y_hat", mean, median, quantile2) %>% 
  mutate(Indices = str_extract_all(variable, "\\d+"), 
         Subject = map_dbl(Indices, \(x) as.integer(x[1])) %>% 
                     factor(levels = levels(data$Subject)),
         Affect = map_chr(Indices, \(x) ifelse(x[2] == 1, "Pos", "Neg")),
         Time_Index = map_dbl(Indices, \(x) as.integer(x[3])))
  
data_predict <- data %>% 
  mutate(DateTime = as_datetime(ymd_hms(paste(as_date(as.double(Day)), 
                                              as.character(Time))))) %>% 
  as_tsibble(key = c(Subject, Affect), 
             index = DateTime) %>% 
  pivot_wider(names_from = Affect, values_from = Score) %>% 
  select(Pos, Neg) %>% 
  drop_na(Pos, Neg) %>% 
  group_by(Subject) %>% 
  mutate(Time_Index = 1:n()) %>% 
  ungroup() %>% 
  pivot_longer(c("Pos", "Neg"), names_to = "Affect", values_to = "Score") %>% 
  left_join(y_hat_summary)

for (i in 1:10) {
  g <- data_predict %>% 
    filter(Subject %in% (10 * (i - 1) + 1):(10 * i)) %>% 
    ggplot(aes(x = DateTime, y = Score)) + 
    geom_line(aes(color = Affect)) + 
    geom_point(aes(color = Affect)) +
    scale_color_manual(values = pos_neg_color) +
    geom_line(aes(y = mean, group = Affect), linetype = "dashed") +
    geom_ribbon(aes(ymin = q5, ymax = q95, group = Affect), alpha = 0.25) +
    geom_hline(yintercept = c(0, 100), color = "forestgreen") +
    scale_y_continuous(limits = c(-20, 120)) +
    scale_x_datetime(breaks = as_datetime(1:7 * 86400),
                     labels = paste("Day", 1:7),
                     limits = as_datetime(c(1, 8) * 86400)) +
    facet_grid(Subject ~ .) 
  
  ggsave(filename = str_glue("plots/mssm2/predict_{from}-{to}.png",
                             from = 10 * (i - 1) + 1, to = 10 * i),
         plot = g, width = 7, height = 14)
}

The fitting recovery results for Subjects 1-20. See the other subject results from this link.

Figure 4: Fitting results for Subjects 1-20. The black dashed lines are the fitted values and the gray areas are 95% CI.

3.2.4 Within- and Between-Subject Reliability

rel_W_draws <- mssm_fit$draws(variables = "rel_W", format = "df") %>%
  select(-.chain, -.iteration, -.draw) %>%
  pivot_longer(cols = everything(), 
               names_to = "variable", values_to = "value") %>% 
  mutate(Indices = str_extract_all(variable, "\\d+"), 
         Subject = map_dbl(Indices, \(x) as.integer(x[1])) %>% 
           factor(levels = levels(data$Subject)),
         Affect = map_chr(Indices, \(x) ifelse(x[2] == 1, "Pos", "Neg")))

g_rel_W <- ggplot(rel_W_draws, aes(x = 1, y = value, fill = Affect)) +
  geom_split_violin() +
  scale_x_continuous(name = NULL, labels = NULL, breaks = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_fill_manual(values = pos_neg_color) +
  facet_wrap(~ Subject, ncol = 10)

ggsave("plots/rel-W.png", g_rel_W, width = 14, height = 14)
Figure 5: The posterior distributions of within subject reliability for PA (red) and NA (blue).
rel_B_draws <- mssm_fit$draws(variables = "rel_B", format = "df") %>%
  mutate(rel_B_diff = `rel_B[1]` - `rel_B[2]`)

g_rel_B <- mcmc_areas(rel_B_draws,
           prob = 0.8,
           prob_outer = 0.99) +
  coord_cartesian(xlim = c(-1.5, 1.5))

ggsave("plots/rel-B.png", g_rel_B)

3.3 Re-parameterize the model

Refer from: https://mc-stan.org/docs/stan-users-guide/multivariate-hierarchical-priors.html

\[ \mathbf{R} = diag.matrix(\boldsymbol{\tau}_R) \times \boldsymbol{\Omega} \times diag.matrix(\boldsymbol{\tau}_R) \]

where \(\boldsymbol{\tau}_R\) is the vector of coefficient scale, \(\boldsymbol{\Omega}\) is the correlation matrix, which can be decomposed by a product of the lower triangle matrix \(\boldsymbol{\Omega_L}\)through Cholesky factorization, \(\boldsymbol{\Omega}_R = \boldsymbol{\Omega_{RL} \Omega_{RL}'}\), for optimization.

\[\begin{align*} \boldsymbol{\tau}_R &\sim Cauchy(loc_{\tau_R}, scale_{\tau_R}) \\ \boldsymbol{\Omega}_{RL} &\sim LKJ.corr.cholesky(\eta_{\Omega_R}) \end{align*}\]

The same setting is applied for \(\mathbf{Q}\).

multilevel-ssm.stan
#include ssm-function.stan

data {
  int<lower=1> N; // number of subjects
  array[N] int<lower=1> T; // number of observation for each subject
  int<lower=1> max_T; // maximum number of observation
  int<lower=1> P; // number of affects
  array[N] matrix[P, max_T] y; // observations 
  array[N] vector[P] m_0; // prior mean of the intial state
  array[N] cov_matrix[P] C_0; // prior covariance of the intial state
}

parameters {
  array[N] vector[P] mu; // ground mean/trane
  array[N] vector[P] theta_0; // initial latent state
  array[N] matrix[P, max_T] theta; // latent states
  array[N] matrix[P, P] Phi; // autoregressive parameters
  
  array[N] cholesky_factor_corr[P] L_Omega_R; 
  array[N] cholesky_factor_corr[P] L_Omega_Q;
  array[N] vector<lower=0>[P] tau_R;
  array[N] vector<lower=0>[P] tau_Q;
  
  vector[P] gamma_mu; // prior mean of the ground mean
  cov_matrix[P] Psi_mu; // prior covariance of the ground mean
  vector[P * P] gamma_Phi; // prior mean of the autoregressive parameters
  cov_matrix[P * P] Psi_Phi; // prior covariance of the autoregressive parameters
  
  vector<lower=0>[P] loc_tau_R; // location of the measurement error
  vector<lower=0>[P] loc_tau_Q; // location of the innovation noise
  vector<lower=0>[P] scale_tau_R; // scale of the measurement error
  vector<lower=0>[P] scale_tau_Q; // scale of the innovation noise
  real<lower=0> eta_Omega_R; // shape of the LKJ prior for the measurement error
  real<lower=0> eta_Omega_Q; // shape of the LKJ prior for the innovation noise
}

transformed parameters{
  array[N] cov_matrix[P] R; // covariance of the measurment error
  array[N] cov_matrix[P] Q; // covariance of the innovation noise
  for (n in 1:N) {
    R[n] = diag_pre_multiply(tau_R[n], L_Omega_R[n]) * diag_pre_multiply(tau_R[n], L_Omega_R[n])';
    Q[n] = diag_pre_multiply(tau_Q[n], L_Omega_Q[n]) * diag_pre_multiply(tau_Q[n], L_Omega_Q[n])';
  }
}

model {
  // level 1 (within subject)
  for (n in 1:N) {
    // when t = 0
    theta_0[n] ~ multi_normal(m_0[n], C_0[n]);
  
    // when t = 1
    theta[n][, 1] ~ multi_normal(Phi[n] * theta_0[n], Q[n]);
    y[n][, 1] ~ multi_normal(mu[n] + theta[n][, 1], R[n]);
    
    // when t = 2, ..., T 
    for (t in 2:T[n]) {
      theta[n][, t] ~ multi_normal(Phi[n] * theta[n][, t - 1], Q[n]);
      y[n][, t] ~ multi_normal(mu[n] + theta[n][, t], R[n]);
    }
  }
  
  
  
  // level 2 (between subject)
  for (n in 1:N) {
    mu[n] ~ multi_normal(gamma_mu, Psi_mu);
    to_vector(Phi[n]) ~ multi_normal(gamma_Phi, Psi_Phi);
    
    
    tau_R[n] ~ cauchy(loc_tau_R, scale_tau_R);
    tau_Q[n] ~ cauchy(loc_tau_Q, scale_tau_Q);
    L_Omega_R[n] ~ lkj_corr_cholesky(eta_Omega_R);
    L_Omega_Q[n] ~ lkj_corr_cholesky(eta_Omega_Q);
    
  }
  
  // the (hyper)priors of parameters are set as the Stan default values
}

generated quantities {
  array[N] matrix[P, max_T] y_hat;
  array[N] matrix[P, P] Tau; 
  array[N] vector[P] rel_W;
  vector[P] rel_B;
  
  for (n in 1:N) {
    // prediction 
    for (t in 1:T[n]) {
      y_hat[n][, t] = mu[n] + theta[n][, t];
    }
    
    // within-subject reliability
    Tau[n] = to_matrix((identity_matrix(P * P) - kronecker_prod(Phi[n], Phi[n])) \ to_vector(Q[n]), P, P);
  
    for (p in 1:P) {
      rel_W[n, p] = Tau[n, p, p] / (Tau[n, p, p] + R[n, p, p]);
    }
  }
  
  // between-subject reliability
  for (p in 1:P) {
    rel_B[p] = Psi_mu[p, p] / (Psi_mu[p, p] + mean(Tau[, p, p]) + loc_tau_R[p]^2);
  }
}
mssm2 <- cmdstan_model("stan/multilevel-ssm2.stan")

data_list <- data %>% 
  #filter(Subject %in% c(1:9, 87)) %>% 
  group_by(Subject) %>% 
  pivot_wider(names_from = Affect, values_from = Score) %>% 
  select(Pos, Neg) %>% 
  drop_na(Pos, Neg) %>% 
  nest() %>% ungroup() %>% 
  mutate(`T` = map_dbl(data, nrow),
         max_T = max(`T`),
         data_padding = pmap(list(data, `T`, max_T), 
                             \(x, y, z) {
                               bind_rows(x, 
                                         tibble(Pos = rep(Inf, z - y),
                                                Neg = rep(Inf, z - y))) %>% 
                                 t()
                             }))


mssm_data <- lst(N = nrow(data_list),
                 `T` = map(data_list$data, nrow),
                 max_T = max(data_list$`T`),
                 P = 2,
                 y = data_list$data_padding, 
                 m_0 = rep(list(c(0, 0)), N), 
                 C_0 = rep(list(diag(c(10^3, 10^3), 2, 2)), N))

mssm2_fit <- mssm2$sample(data = mssm_data, 
                          seed = 20240221,
                          chains = 8,
                          parallel_chains = 8,
                          iter_sampling = 1500,
                          refresh = 1000, 
                          show_messages = FALSE)
mssm2_fit$save_object(file = "stan/multilevel-ssm2-fit.RDS")
saveRDS(mssm2_fit$summary(), "stan/multilevel-ssm2-summary.RDS")
Data summary
Name data.frame(Rhat = as.doub…
Number of rows 31661
Number of columns 2
_______________________
Column type frequency:
numeric 2
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Rhat 3002 0.91 1.55 0.46 1.00 1.31 1.4 1.56 5.12 ▇▁▁▁▁
ESS 3002 0.91 31.39 172.43 8.37 14.76 19.0 25.70 12032.13 ▇▁▁▁▁
rel_W_draws_2 <- mssm2_fit$draws(variables = "rel_W", format = "df") %>%
  select(-.chain, -.iteration, -.draw) %>%
  pivot_longer(cols = everything(), 
               names_to = "variable", values_to = "value") %>% 
  mutate(Indices = str_extract_all(variable, "\\d+"), 
         Subject = map_dbl(Indices, \(x) as.integer(x[1])) %>% 
           factor(levels = levels(data$Subject)),
         Affect = map_chr(Indices, \(x) ifelse(x[2] == 1, "Pos", "Neg")))

g_rel_W_2 <- ggplot(rel_W_draws_2, aes(x = 1, y = value, fill = Affect)) +
  geom_split_violin() +
  scale_x_continuous(name = NULL, labels = NULL, breaks = NULL) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_fill_manual(values = pos_neg_color) +
  facet_wrap(~ Subject, ncol = 10)

ggsave("plots/rel-W-2.png", g_rel_W_2, width = 14, height = 14)

rel_B_draws_2 <- mssm2_fit$draws(variables = "rel_B", format = "df") %>%
  mutate(rel_B_diff = `rel_B[1]` - `rel_B[2]`)

g_rel_B_2 <- mcmc_areas(rel_B_draws_2,
           prob = 0.8,
           prob_outer = 0.99) +
  coord_cartesian(xlim = c(-1.5, 1.5))

ggsave("plots/rel-B-2.png", g_rel_B_2)
Figure 6: The posterior distributions of within subject reliability for PA (red) and NA (blue).

The posterior distributions of within subject reliability for PA (red) and NA (blue).

However, the results are worse than the original one. The main reason is the MCMC samples had not converged yet.

References

Koval, P., Pe, M. L., Meers, K., & Kuppens, P. (2013). Affect dynamics in relation to depressive symptoms: Variable, unstable or inert? Emotion, 13(6), 1132–1141. https://doi.org/10.1037/a0033579
Molenberghs, G., Laenen, A., & Vangeneugden, T. (2007). Estimating reliability and generalizability from hierarchical biomedical data. Journal of Biopharmaceutical Statistics, 17(4), 595–627. https://doi.org/10.1080/10543400701329448
Schuurman, N. K., & Hamaker, E. L. (2019). Measurement error and person-specific reliability in multilevel autoregressive modeling. Psychological Methods, 24(1), 70–91. https://doi.org/10.1037/met0000188
Vangeneugden, T., Laenen, A., Geys, H., Renard, D., & Molenberghs, G. (2005). Applying concepts of generalizability theory on clinical trial data to investigate sources of variation and their impact on reliability. Biometrics, 61(1), 295–304. https://doi.org/10.1111/j.0006-341X.2005.031040.x